home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / sum / test.c < prev    next >
Encoding:
C/C++ Source or Header  |  2001-11-01  |  5.0 KB  |  235 lines
  1. /* sum/test.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <stdlib.h>
  24. #include <stdio.h>
  25. #include <gsl/gsl_math.h>
  26. #include <gsl/gsl_test.h>
  27. #include <gsl/gsl_sum.h>
  28.  
  29. #include <gsl/gsl_ieee_utils.h>
  30.  
  31. #define N 50
  32.  
  33. void check_trunc (double * t, double expected, const char * desc);
  34. void check_full (double * t, double expected, const char * desc);
  35.  
  36. int
  37. main (void)
  38. {
  39.   gsl_ieee_env_setup ();
  40.  
  41.   {
  42.     double t[N];
  43.     int n;
  44.  
  45.     const double zeta_2 = M_PI * M_PI / 6.0;
  46.  
  47.     /* terms for zeta(2) */
  48.  
  49.     for (n = 0; n < N; n++)
  50.       {
  51.     double np1 = n + 1.0;
  52.     t[n] = 1.0 / (np1 * np1);
  53.       }
  54.  
  55.     check_trunc (t, zeta_2, "zeta(2)");
  56.     check_full (t, zeta_2, "zeta(2)");
  57.   }
  58.  
  59.   {
  60.     double t[N];
  61.     double x, y;
  62.     int n;
  63.  
  64.     /* terms for exp(10.0) */
  65.     x = 10.0;
  66.     y = exp(x);
  67.  
  68.     t[0] = 1.0;
  69.     for (n = 1; n < N; n++)
  70.       {
  71.     t[n] = t[n - 1] * (x / n);
  72.       }
  73.  
  74.     check_trunc (t, y, "exp(10)");
  75.     check_full (t, y, "exp(10)");
  76.   }
  77.  
  78.   {
  79.     double t[N];
  80.     double x, y;
  81.     int n;
  82.  
  83.     /* terms for exp(-10.0) */
  84.     x = -10.0;
  85.     y = exp(x);
  86.  
  87.     t[0] = 1.0;
  88.     for (n = 1; n < N; n++)
  89.       {
  90.     t[n] = t[n - 1] * (x / n);
  91.       }
  92.  
  93.     check_trunc (t, y, "exp(-10)");
  94.     check_full (t, y, "exp(-10)");
  95.   }
  96.  
  97.   {
  98.     double t[N];
  99.     double x, y;
  100.     int n;
  101.  
  102.     /* terms for -log(1-x) */
  103.     x = 0.5;
  104.     y = -log(1-x);
  105.     t[0] = x;
  106.     for (n = 1; n < N; n++)
  107.       {
  108.     t[n] = t[n - 1] * (x * n) / (n + 1.0);
  109.       }
  110.  
  111.     check_trunc (t, y, "-log(1/2)");
  112.     check_full (t, y, "-log(1/2)");
  113.   }
  114.  
  115.   {
  116.     double t[N];
  117.     double x, y;
  118.     int n;
  119.  
  120.     /* terms for -log(1-x) */
  121.     x = -1.0;
  122.     y = -log(1-x);
  123.     t[0] = x;
  124.     for (n = 1; n < N; n++)
  125.       {
  126.     t[n] = t[n - 1] * (x * n) / (n + 1.0);
  127.       }
  128.  
  129.     check_trunc (t, y, "-log(2)");
  130.     check_full (t, y, "-log(2)");
  131.   }
  132.  
  133.   {
  134.     double t[N];
  135.     int n;
  136.  
  137.     double result = 0.192594048773;
  138.  
  139.     /* terms for an alternating asymptotic series */
  140.  
  141.     t[0] = 3.0 / (M_PI * M_PI);
  142.  
  143.     for (n = 1; n < N; n++)
  144.       {
  145.     t[n] = -t[n - 1] * (4.0 * (n + 1.0) - 1.0) / (M_PI * M_PI);
  146.       }
  147.  
  148.     check_trunc (t, result, "asymptotic series");
  149.     check_full (t, result, "asymptotic series");
  150.   }
  151.  
  152.   {
  153.     double t[N];
  154.     int n;
  155.  
  156.     /* Euler's gamma from GNU Calc (precision = 32) */
  157.  
  158.     double result = 0.5772156649015328606065120900824; 
  159.  
  160.     /* terms for Euler's gamma */
  161.  
  162.     t[0] = 1.0;
  163.  
  164.     for (n = 1; n < N; n++)
  165.       {
  166.     t[n] = 1/(n+1.0) + log(n/(n+1.0));
  167.       }
  168.  
  169.     check_trunc (t, result, "Euler's constant");
  170.     check_full (t, result, "Euler's constant");
  171.   }
  172.  
  173.   {
  174.     double t[N];
  175.     int n;
  176.  
  177.     /* eta(1/2) = sum_{k=1}^{\infty} (-1)^(k+1) / sqrt(k)
  178.  
  179.        From Levin, Intern. J. Computer Math. B3:371--388, 1973.
  180.  
  181.        I=(1-sqrt(2))zeta(1/2)
  182.         =(2/sqrt(pi))*integ(1/(exp(x^2)+1),x,0,inf) */
  183.  
  184.     double result = 0.6048986434216305;  /* approx */
  185.  
  186.     /* terms for eta(1/2) */
  187.  
  188.     for (n = 0; n < N; n++)
  189.       {
  190.     t[n] = (n%2 ? -1 : 1) * 1.0 /sqrt(n + 1.0);
  191.       }
  192.  
  193.     check_trunc (t, result, "eta(1/2)");
  194.     check_full (t, result, "eta(1/2)");
  195.   }
  196.  
  197.   exit (gsl_test_summary ());
  198. }
  199.  
  200. void
  201. check_trunc (double * t, double expected, const char * desc)
  202. {
  203.   double sum_accel, prec;
  204.  
  205.   gsl_sum_levin_utrunc_workspace * w = gsl_sum_levin_utrunc_alloc (N);
  206.   
  207.   gsl_sum_levin_utrunc_accel (t, N, w, &sum_accel, &prec);
  208.   gsl_test_rel (sum_accel, expected, 1e-8, "trunc result, %s", desc);
  209.  
  210.   /* No need to check precision for truncated result since this is not
  211.      a meaningful number */
  212.  
  213.   gsl_sum_levin_utrunc_free (w);
  214. }
  215.  
  216. void
  217. check_full (double * t, double expected, const char * desc)
  218. {
  219.   double sum_accel, err_est, sd_actual, sd_est;
  220.   
  221.   gsl_sum_levin_u_workspace * w = gsl_sum_levin_u_alloc (N);
  222.  
  223.   gsl_sum_levin_u_accel (t, N, w, &sum_accel, &err_est);
  224.   gsl_test_rel (sum_accel, expected, 1e-8, "full result, %s", desc);
  225.   
  226.   sd_est = -log10 (err_est/fabs(sum_accel));
  227.   sd_actual = -log10 (DBL_EPSILON + fabs ((sum_accel - expected)/expected));
  228.  
  229.   /* Allow one digit of slop */
  230.  
  231.   gsl_test (sd_est > sd_actual + 1.0, "full significant digits, %s (%g vs %g)", desc, sd_est, sd_actual);
  232.  
  233.   gsl_sum_levin_u_free (w);
  234. }
  235.